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Abstract. We investigate the distribution of work performed on a Brownian particle 
in a time-dependent asymmetric potential well. The potential has a harmonic 
component with time-dependent force constant and a time-independent logarithmic 
barrier at the origin. For arbitrary driving protocol, the problem of solving the Fokker- 
Planck equation for the joint probability density of work and particle position is reduced 
to the solution of the Riccati differential equation. For a particular choice of the driving 
protocol, an exact solution of the Riccati equation is presented. Asymptotic analysis of 
the resulting expression yields the tail behavior of the work distribution for small and 
large work values. In the limit of vanishing logarithmic barrier, the work distribution 
for the breathing parabola model is obtained. 
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1. Introduction 

Stochastic thermodynamics is an advancing field with many applications to small 
systems of current interest [H l2l [3] . Work performed on a small system by an external 
driving becomes a stochastic variable because of the strong infiuence of fiuctuations 
mediated by the environment. In a standard setting, a particle moves in a thermal 
environment and experiences a time-dependent external force. The particle position by 
itself is a stochastic process, say X(t). Any single trajectory of the particle in a time 
interval [0,t] yields a single value of the work W(t) done on the particle. The work 
W(t) thus is a functional of the position process X(t'), < t' < t, and is distributed 
with a probability density function (PDF) p{w;t). The probability p{w;t)dw that 
the work W(t) falls into an infinitesimal interval {w, w + dw) equals the probabilistic 
weight of all trajectories giving work values in that interval. An important aspect 
of stochastic thermodynamics is that relations between work distributions for forward 
and reversed driving protocols as well as averages over functions of the work allow one 
to determine thermodynamic potential differences between equilibrium states of small 
systems. Perhaps most widely known is the Jarzynski equality [4], which relates the free 
energy difference between two equilibrium configurations of a system to the average 



where = k-^T is the thermal energy. In computer simulations and experiments, 
sampled values of w lie typically within one or two standard deviations of the maximum 
of p{w]t), while values most important in Eq. ([1]) are those near the maximum of 
exp{—l3w)p{w]t). The corresponding regimes overlap significantly only if the function 
exp(— does not change much over one standard deviation. If this is not the case, the 
relevant contributions to the integral come from those rare trajectories with work value 
belonging to the tails of p{w; t). In experiments these rare trajectories are almost never 
observed and even in simulations it is difficult to generate them with sufficient statistical 
weight. Hence, to evaluate averages as in Eq. ([T]), information on the tails of the PDF 
is essential. An important part of this study is to gain insight into the asymptotic PDF 
behavior of the work performed on a Brownian particle in an asymmetric time-dependent 
potential well. 

We will calculate the characteristic function for the work in a simple setting which, 
however, may be realized in experiments [HlEllTj. In this setting, an overdamped motion 
of a Brownian particle is considered in the logarithmic-harmonic potential 



where the parameter (7 > specifies the strength of the logarithmic repulsive wall near 
the origin, and k{t) is a time- dependent force constant. In the deterministic limit, 
i.e., in absence of thermal noice, the particle moves along the positive x-axis as driven 
by the time- dependent force —dU{x,t)/dx (without inertia). Taking into account the 




(1) 




g > 0, X > 0. 



(2) 



Work distribution in logarithmic-harmonic potential 



3 



thermal noise, the combined process {X(t), W(t)} is described by the system of Langevin 
equations 



where D quantifies the strength of the noise and B(t) is the standard Wiener process. 
Specifically, in the case of a thermal environment, the noise strength is proportional to 
the temperature, D = k-oT (when the particle mobility is set to one). 

While the work process has not yet been studied for the logarithmic-harmonic 
potential, an exact solution of the Fokker-Planck equation for the position process was 
obtained |8]. In the following, we recover this solution from our Lie algebraic approach. 
The solvability of this problem stems from the fact that the operators entering the 
Fokker-Planck equation form a Lie algebra [HI ID]- If one considers the Fokker-Planck 
equation for the joint PDF of position and work, the corresponding differential operators 
no longer form a Lie algebra. However, if one starts with the Fokker-Planck equation for 
the joint PDF and performs a Laplace transformation with respect to the work variable 
w, a Lie algebra is obtained. Solution of the Fokker-Planck equation then provides 
the characteristic function for the work process and the tails of the work PDF can be 
extracted using asymptotic analysis of Laplace transforms. 

The work PDF for the problem of the "breathing parabola" [g = 0) has been studied 
analytically in Refs. [11], [12], [13], and [H]. In [11], [12], the authors considered an 
expansion around a single trajectory attributed to a prescribed rare value of the work 
and derived asymptotic results for the tails of the work PDF in the small temperature 
limit. The solution reported in [13] is formally exact for arbitrary protocol k{t). Explicit 
results are given in the limit of slow driving, where the process is close to a quasi-static 
equilibrium and the work PDF can be approximated by a Gaussian. In Ref. [Hj the 
work-weighted propagator was derived by the path integral method. Another closely 
related setting, where the work PDF can be calculated analytically, is for a parabolic 
potential with a time- dependent position of the minimum ( "sliding parabola" ) [H] , [I5] , 
[T6] . [T7] . [T8] . The present work broadens the list of few exact results in this field. 

2. Solution of Fokker-Planck equation for arbitrary protocol 

2.1. Green's function for logarithmic potential 

An important auxiliary quantity in deriving all subsequent results is the Green's function 




(3) 




(4) 




6{x — Xo) , Xo > , 



(5) 



which represents the solution of the Fokker-Planck equation 




(6) 



Work distribution in logarithmic-harmonic potential 



4 



for the diffusion in the time- independent logarithmic potential with initial condition 
q{x] 0|xo) = 6{x — Xq). The explicit form of the solution is [HI 120] 

expl 



q{x;t\xo) 



Xo f x_ 
2t \xo 



X2 I ^2 



At 



) 'm . 



where I^i-) is the modified Bessel function of order u, and 



2 



(7) 



(8) 



measures the strength of the logarithmic potential in relation to the intensity of the 
thermal noise. In all subsequent results the parameter g enters solely through u defined 
in Eq. ([8]). 

Equation ([7]) is the unique norm-preserving solution of the diffusion problem in the 
domain a; > 0, i.e., the probability current at a; = vanishes. Therefore, performing the 
limit g ^ and using /_i/2(-2) = cosh{z) / ^/z, we obtain the standard solution 

for free diffusion with a reflecting boundary at a; = [2T] . 



2.2. Joint Green's function for work and position 

Let us denote by p{x,w;t\xo,0) the joint PDF for the process {X(t),W(t)} given that 
at time t = the particle is at position xq, xq > 0, and no work has been done yet, 

p{x,W]0\Xo,0) = 6{x — Xq)6{w) , Xq > . (9) 

The time evolution of the joint PDF is given by the Fokker-Planck equation 



d 

—p{x,w;t\xo,0) 
at 



D— -^{l- k{t)x] - -k{t)x' — 
dx^ dx \x / 2 dw 



p{x, w; t\xo, 0) . 



(10) 



The differential operators on the right-hand side do not exhibit closed commutation 
relations. However, after performing the two-sided Laplace transformation 



p{x,C.;t\xo) 



+ 00 

dwe" 

-oo 



"p{x, w; t\xo, 0) 



the Fokker-Planck equation (ITU]) assumes the form 

|-p(x, t|xo, 0) = \dJo + 2k{t)Ji - A^k{t)J2 + + I) k{t) p{x, ^; t|xo, 0) 
at I 

where the differential operators 

_ g d 1 
D dx X ^ 



Jn 



^ If d 



J2 



-X 



satisfy the closed commutation relations 

[^07 ^l] = ^0 ) [^05 = Jl 



:iii 



(12) 



(13) 



(14) 



This allows us to apply the Lie algebraic method, as discussed, e.g., in Refs. [9l 123] . 
First, we write the solution of Eq. (fT2|) in the factorized form 



p{x,i]t\xo) 



exp 



[v + l) dt'k(t') 
Jo 



exp( 62(^)^2) exp(6i(t)Ji) exp[bo{t)Jo)S{x-Xo) (15) 
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where the time-dependent coefficient 62 (^) is obtained by solving the Riccati differential 
equation 

Ht) = ^hl{t) + 2k{t)h2{t) - Aik{t) , 62(0) = . (16) 
Knowing 62 (^) the other coefficients are given by 

61 (t) = 2 f dt'k{t') + D [ dt'h^it'), bo{t) = D [ dt'exp[bi{t')] . (17) 
Jo Jo Jo 

In the last step, we evaluate, using Eqs. ([5]) and ([7]), the action of the operator 

exp{bo{t)Jo) on the delta function in Eq. (|T5]) . and subsequently apply to the 

corresponding result the two remaining exponential operators in Eq. f|T5|) . This yields 



p{x,C]t\xo) =exp(^jdt'k{t')-^uDjdt'b2{t') + ^b2{t)x^^q{xe^^ . (18) 

In the derivation we have utilized the operator identity 

exp (?7x^ ) /(x) = /(2;exp(?7)) . (19) 



dx ^ 

The exact solution (fTSl) is the central result of the present section and constitutes the 
starting point of all subsequent analyses. 



3. PDF of particle position and its long-time asymptotics 



After integrating the joint PDF p{x,w]t\xQ,0) over the work variable, the transition 
PDF p{x;t\xo) for the particle coordinate is obtained. Equivalently, the w-integration 
is accomplished by evaluating the result (fT8|) at ^ = [cf. Eq. (fTT!) ]. Notice that the 
variable ^ enters the solution (fT8|) only through the Riccati equation (fT6!l . When taking 
^ = 0, this equation reduces to the Bernoulli differential equation, where the unique 
solution satisfying 62(0) = is the trivial 01 
Eq. ( |T8l) are then given by 



b,{t) = 2 dt'k{t'), bo{t) 
Jo 



D 



62 (t) = 


The remainin 






/ dt' exp 


2 fdt"k{t") 


Jo 


Jo 



Hence the PDF for the particle position read 



p{x]t\xQ) 



X 

This finding is valid for an arbitrary driving protocol k{t). If k{t) is a positive constant, 
say k{t) = fco > 0, then the system approaches the Gibbs canonical equilibrium at long 
times. If the constant force constant is superimposed with a periodically oscillating 
component, a gradual formation of a nontrivial steady state occurs. In this steady 

I This result agrees with Eq. (19) in Ref. [24 , where it has been derived in connection with a diffusion 
problem with logarithmic factors in drift and diffusion coefhcients. 



3,2gfei(t) 



+ 4 



xxq e 2 



lit) 



2bo{t) 



(20) 



(21) 
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state, the PDF does not depend on the initial condition xq and, for any given x > is 
a periodic function of time with the fundamental period given by that of k{t). 
To exemplify the PDF in the steady state, let us take 

k{t) = ko + ki sm{ut) , ko>0. (22) 

The asymptotic analysis of Eq. ( 12T|) for long times requires the evaluation of the limit 

hm— -, a<l. (23) 

t^oo bo[t) 

If a < 1, the limit exists and, using L'Hopital's rule, it equals zero. Hence for any 
finite x and xq, the argument z = xxqc^^^^^^"^ / {2bo{t)) in the Bessel function appearing 
in Eq. fl^T]) becomes small for large t and we can write Iu{z) ~ {^zY /Tlu + 1). If a = 1, 
the limit does not exist and e'''^^^^ /bo{t) ~ 1/ f{t) for t oo^ where 

f{t) = D exp cos(a;t) /„, ^, ^. , fco > . (24) 



Accordingly, for t — oo 



p(x;tM~j>as(x;t) = ^;^^(^-^J (2 J 



X\2 1 



.2/ m\ ■ 

In the limit /ci — > or a; — )■ 0, f{t) — )■ D/{2ko), and p(x; t|a;o) approaches the Gibbs 
equilibrium distribution. 

Finally, the limit (7 — )■ in Eq. (pTj) [Eq. (!25|) ] yields the exact transition PDF (exact 
time-asymptotic PDF) for the breathing parabola model with reflecting boundary at the 
origin. Since the parameter g enters only via z/ in Eq. ([8]), (7 — ?■ corresponds to — t- — ^ 
in Eqs. (12T]) and (|25|). 



4. Work fluctuations 



. i . Characteristic functions 

By integration of the joint PDF in Eq. f lTSj) over the spatial variable x, we obtain the 
characteristic function for the work done on the particle during the time interval [0,t]. 

Let us first consider the particle dynamics conditioned on the initial position xq. 
In this case the characteristic function for the work reads 



dxp{x,^; t\xo) 



Carrying out the integration we find 



^^■,t\xo) 



2ef'i(*)-|^/o'it'fc2(t') 
2e''iW -6o(t)62(t) 



u+l 



exp 



2 J 2e''iW -6o(t)&2(t) . 



(26) 



(27) 



A physically more important situation is, when the particle coordinate is initially 
equilibrated with respect to the initial value fc(0) = fco > of the force constant. In 
order to obtain the characteristic function for this situation we have to integrate over 
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Xq the product of the characteristic function $(^; t\xo) and the equihbrium PDF pas for 
f{t) = D/{2ko) given in Eq. ([25]). The result is 

Note that for C,= (3, the average in Eq. ([T]) equals ^{(3; t). Equation ( l28l) is valid for an 
arbitrary driving protocol k{t). The Laplace variable ^ enters $(^; t) implicitly through 
the function &2(^) [cf- the Riccati equation f|T6|) ] and through the functions &i(t) and 
6o(t) [cf. Eq. (HZ])]. 

In the limit (? — (z/ — > — |), Eqs. (I27II28P give the corresponding characteristic 
functions for the breathing parabola model with reflecting boundary at x = 0. These 
characteristic functions are also valid for the breathing parabola model without reflecting 
boundary, if rather obvious changes are made of the meaning of the initial coordinate 
xo in Eq. ( 1271) . and of the initial Gibbs equilibrium state underlying Eq. ( 128]) . In 
the breathing parabola model without reflecting boundary, Eq. (1271) is valid for xq G 
(— oo, +oo) and Eq. (128]) corresponds to the initial Gibbs equilibrium in the parabolic 
potential U{xo) = kQx'^/2. The equivalence of the characteristic functions for the 
problems with and without reflecting boundary is due to the symmetry of the parabolic 
potential, which implies that the work done on the particle that crosses the origin is the 
same as the work done on the particle reflected at the origin. This reasoning can be 
supported by an independent calculation if one notices that both models, the present 
model with logarithmic-harmonic potential and the breathing parabola one, have the 
same operator algebra. 

4-2. Simple example 

A driving protocol, where an explicit solution of the Riccati equation (1161) in terms of 
elementary functions can be given, is 

1 + 7t 

Notice that the same protocol was considered before in Refs. [HI [T2]. According to 
Eq. fl29]) . the potential well widens with time and hence the work done on the particle 
is negative for any x > 0. The solution of the Riccati equation flTB]) reads 

2 d 



62(t) = -;^^ln^l+7t)^'^''°^^^-^^^^l 



^1 + ^t)^iO _ (2^0 + 7) + A{0 



(30) 



(2fco + 7)-A(0. 
where 

A{i) = ^(2ko + jy-8kojD^ . (31) 

For simplicity, we take 7 = 1 and /cq = 1 in the following. After calculating boit) and 
bi{t) from Eqs. (JT7]) . the characteristic function in Eq. (128]) becomes 

( 2,4(0(1 + t)il'^-^<«l y" 

I 1(1 + i)-'"^' + 1] + (3 - 2D() [(1 + t)-MO - 1] I ■ I 
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For ^ = 1/D in particular, we obtain t) = (1 + t)^"*"^ which exemphfies the 

Jarzynski equahty for the driving protocol fl29|) . 

Successive derivatives of the characteristic function with respect to ^ evaluated at 
^ = yield the cumulants of the work distribution. The mean work done on the particle 
during the time interval [0, t] is 



61n(l+t) + 



+ 3t^ + 3t 



(33) 



It is a monotonically decreasing function of t, where for small times, the decrease is 
linear, while in the long-time limit, it is logarithmic. For the variance we find 

(W^(t))-(W(t))^ = (z. + l)(^^y '' + y+/^ [t^ + 3t^ + 3t + 24(l+t)3ln(l+t)] .(34) 

This increases monotonically, where the increase is quadratic for small times and 
logarithmic for long times. The strength g of the logarithmic potential barrier enters 
the above formulae only through the multiplicative prefactor {u + 1). This holds true 
for all cumulants. For stronger repulsion, the particle predominantly diffuses in a region 
further away from the origin. The decrease of its typical potential energy results in a 
larger absolute value of the mean work. At the same time, the width of the work PDF 
increases, since the initial particle position is sampled from a broader Gibbs distribution. 

Equation ( 132|) entails the complete information about the work distribution p{w; t). 
In particular it allows one to derive the tails of the work PDF for both w — )■ 0~ and 
w — )■ — oo without carrying out the inverse Laplace transformation of $(^;t). 

The asymptotics of p{w;t) for — )■ 0~ is related to the asymptotics of $(^;t) for 
^ — 7- — oo, which follows from (!32l) . 

By taking the inverse Laplace transform of this asymptotic form (cf. Ref. [25]), we obtain 
the parabolic cylinder function (cf. Ref. [26]) with argument proportional to \J D /\w\. 
Considering the limit of large arguments of this function [26] , we find 



p{w; t) ~ ci(t) ( ^ ) e'^^^^'i^ , w ^ 0" , (36) 



where 

- 2 I — / , . 3 



e,W=^_,„(l + <)j , e.W = -^-(l+«).(^l^;^i|^j . (37) 

For any g and any t, the PDF almost vanishes in an interval {wdt), 0), where its width 
|tfc(t)| is controlled by the "damping constant" C2{t). The width increases both with 
time and the strength of the logarithmic potential. This can be understood from the 
fact that any trajectory yielding a small (absolute) value of the work, must necessarily 
depart from a position close to the origin and remain in its vicinity during the whole 
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2 4 6 8 10 5 10 15 20 

t t 

Figure 1. Positions of tlie singularity in the expression (|38| (left panel), and the 
prefactor ((39)) (right panel) as functions of time. These functions control, through 
Eq. pO|) . the large \w\ asymptotics of the work PDF. In the right panel we have taken 
D = 1 and g = I. 



time interval [0,t]. The probabilistic weight of such trajectories decreases with both t 
and g. 

The asymptotics of the work PDF p{w; t) for w — oo (at fixed t) is determined 
by the expansion of the characteristic function at that ^o{t)y which gives the 

singularity of lying closest to its domain of analyticity [22]. To find ^o(^), we 

numerically solved the transcendental equation l/<l>(^o(^);^) = 0. In the vicinity of the 
singularity, 

m t) ~ r(t) (-^^——] ^ , e ^ Ut) , (38) 



with 



r{t) = {-Dr\ hm (e - ut)r' m t) . m 



From this result we obtain [22] 

P<"'''~^r(7?T)(lf)''''"'''*''^ • <«' 

where D^o{t) is a real, positive, decreasing function of t, cf. Fig. [H For small t, the work 
PDF is very narrow {DC,o{t) large). With increasing t, the weight of the trajectories 
yielding large (absolute) values of the work increases. This is reflected by the decrease 
of -D,^o(^)lil Contrary to the function C2(t) in ( l36l) . D^o{t) does not depend on the 
strength of the logarithmic potential. The parameter g enters only the pre-exponential 
factor in Eq. ( l40l) . 

§ For t = 2 and g = we obtain DCo(2)^1.827, r(2)/r(l/2)=1.021, which is in perfect agreement with 
Eq. (5.30) in [12]. 
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w 



Figure 2. a) Simulated work PDF in comparison with the asymptotic behavior 
predicted by Eq. {\w\ large, sohd line) and Eq. (HT|) {\w\ small, dashed line) 

for parameters g = 1.5, D = 1 and t = 1. In the simulations 10^ trajectories were 
generated with a time step At — 0.001 (adapted when the particle is near the origin, see 
text), b) Semi-logarithmic plot of simulated p(w;t) vs. w (circles), demonstrating the 
agreement with Eq. (solid line) for large \w\. c) Semi- logarithmic plot of simulated 
p{w; t) vs. w (circles) in comparison with the first leading term of the asymptotic 
expansion for small \w\ (Eq. ((36)) . solid line), and when including the second leading 
term according to Eq. ([4T|) (dashed line). 



In order to verify the exact asymptotic expansions of the work PDF, we have 
performed extensive Langevin dynamics simulations using the Heun algorithm [27J for 
several sets of parameters and different time intervals. A typical PDF together with 
the predictions for its asymptotic behavior according to Eqs. ( |36l) and (HOl) is shown 
in Fig. [21 In order to avoid nonphysical negative values of the particle position in the 
numerics, which can originate from a fixed time discretization, we have implemented a 
time-adapted Heun scheme. If a negative (attempted) coordinate along a trajectory is 
generated, the time step At is reduced until the attempted particle position is positive. 
To allow for a better comparison of the analytical findings with the simulated data in 
Fig. [2] for small \w\, we have derived also the second leading term in the asymptotic 
expansion for \w\ — )■ 0. After somewhat lengthy but straightforward calculation, we 
obtain 



\w\ 

15 



where Ci{t) and C2(t) are given in Eq. ( ISTjl . 




4 



t) 



\w\ 

15 



(41) 
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5. Concluding remarks 

Based on a Lie algebraic approach we succeeded to derive Eq. ( fT8l) for the joint PDF of 
work and position for a Brownian particle in a time-dependent logarithmic-harmonic 
potential. In order to derive explicit results from Eq. ( JTSl) for a given protocol, 
the Riccati equation ( IT6l) needs to be solved. This nonlinear differential equation is 
equivalent to the linear second-order differential equation [28] 

y{t)-^Ht)y{t)-^D^Ht)y{t) = , ^(0) = 0. (42) 

Specifically, if y{t) solves fH2l) . then the logarithmic derivative 

^'^'^ = -dW) ' ^ ^ 

is the solution of Eq. (fTB]) . Hence the characteristic function f l28p can be expressed in 
terms of the function y{t). 

The solution of f H2|) for several reasonable driving protocols, e.g., for k{t) = 
fcoexp(±7t), k{t) = fco + kit, or k{t) = kit"-, can be written in terms of higher 
transcendental functions. Corresponding results are quite involved and will be published 
elsewhere. Here we have focused on the simple protocol (1291) which should exemplify 
typical asymptotic features of the work PDF for monotonia driving. Notice that, if 
k{t) > and ^ — )■ cxd along the real axis, one can use the WKB approximation and derive 
a generic expression for y{t) valid for any protocol, i.e., also a generic approximative 
expression for the work characteristic function (128|1 . 

For non-monotonic driving protocols the work can assume any real value. Then 
the work PDF has the support (— oo, +oo) and its two-sided Laplace transform will be 
analytic within a stripe parallel to the imaginary axis. The w — )■ -|-oo {w — )■ — oo) tail 
of the work PDF is determined by the singularity, which is closest to the stripe on its 
left (right) side. We hence expect an asymptotics 



Piw;t)r^- ^^^'^^ , (M]\-Di±it)^4 , w^±oo, (44) 



where the coefficients ^±{t), r±{t), depend on the driving protocol k{t). Periodic driving 
protocols play an important role in the analysis of Brownian motors. A deeper analysis 
of the work PDF for this class of protocols seems to be worthy for further study. 
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